Introduction

Nous souhaitons réalisé l’étude d’une série temporelle et faire des prévisions sur celle-ci.

Cette série temporelle est le trafic mensuel d’une Compagnie aérienne de janvier 2011 à août 2019.

Nos prévisions portent sur les 8 mois de l’année 2019

Représentation graphique de la série.

Import des données

Import de la base, on sélectionne la colonne des valeurs de trafic.

library(readr)
data <- read_delim("Trafic-voyageurs.csv", 
    delim = ";", locale = locale(encoding = "ISO-8859-1"))
data_value <- data[,2]
summary(data)
##     dates               trafic      
##  Length:104         Min.   :220876  
##  Class :character   1st Qu.:297154  
##  Mode  :character   Median :355178  
##                     Mean   :354651  
##                     3rd Qu.:407331  
##                     Max.   :505190

Création de notre série temporelle

Nous représentons nos données sous forme de série temporelle.

Une série temporelle est un ensemble de métrique mesurée sur des intervalles de temps réguliers.

Création de la série chronologique avec la librairie TSstudio :

  • start : année de début,
  • frequency : La fréquence est mensuelle donc on sélectionne une fréquence de 12/an
library(TSstudio)
data_ts <- ts(data_value, start = 2011, frequency = 12)
plot_1_TimeSeries(data_ts)
## 
## Attachement du package : 'plotly'
## L'objet suivant est masqué depuis 'package:ggplot2':
## 
##     last_plot
## L'objet suivant est masqué depuis 'package:stats':
## 
##     filter
## L'objet suivant est masqué depuis 'package:graphics':
## 
##     layout

Séparation jeux de données

data_ts_train <-
  window(data_ts, start = c(2011, 1), end = c(2018, 12))
data_ts_test <- window(data_ts, start = c(2019, 1), end = c(2019, 8))

names(data)[1] <- "ds"
names(data)[2] <- "y"
data_train <- data[1:96, ]
data_test <- data[97:104, ]


plot(data_ts, xlim = c(2011, 2020))
lines(data_ts_test, col = 3)
legend(
  "topleft",
  lty = 1,
  col = c(1, 3),
  legend = c("Série chronologique Train", "Série chronologique Test")
)

On observe une forte tendance, et un paterne semble se répété, s’agit-il d’une saisonnalité ?

Analyse de notre série temporelle

Pour commencer nous allons nous intéresser à savoir si notre série est stationnaire ou non.

La stationnarité d’une série signifie que le processus qui génère la série ne change pas dans le temps. Cela ne veut pas dire que la série ne change pas dans le temps, mais que la façon dont elle change, n’est pas modifié dans le temps.

Testons si la série est stationnaire :

Par la suite la library tsutils permet de vérifier automatiquement s’il y a une tendance et une saisonnalité grâce au test effectué par la fonction seasplot().

library(tseries)
adf.test(data_ts) #p-value <0.5 => on ne rejete pas H0 => non stationnaire
## Warning in adf.test(data_ts): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_ts
## Dickey-Fuller = -5.4857, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(data_ts)
## Warning in kpss.test(data_ts): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  data_ts
## KPSS Level = 2.0803, Truncation lag parameter = 4, p-value = 0.01
library(tsutils)
## Registered S3 method overwritten by 'greybox':
##   method     from
##   print.pcor lava
seasplot(data_ts) #TRUE | TRUE

## Results of statistical testing
## Evidence of trend: TRUE  (pval: 0)
## Evidence of seasonality: TRUE  (pval: 0)

Donc notre série est bien non-stationnaire.

Décomposition de la série chronologique :

Généralement les séries chronologiques sont divisées en plusieurs composantes qui peuvent être additionnées ou multipliées entre elles.

Additive : Ts = Tendance + Saisonnalité+ Erreur

Multiplicative : Ts = Tendance * Saisonnalité * Erreur

Donc chaque point de notre série temporelle peut être exprimer par :

  • Tendance \(Tt\) : qui correspond à l’indicateur de si nos données changent sur le long terme ou non, elle permet de savoir si la série se déplace de manière fluide ou progressivement.

  • Saisonnalité \(St\) : décrit un pattern qui se reproduit sur un intervalle de temps régulier (ici une année), c’est à dire que chaque année, le même pattern est présent.

  • Erreur \(ϵt\) : ou résidus correspond aux “restes” quand on a supprimé la tendance et la saisonnalité. Ces valeurs peuvent s’expliquer par des évènements aléatoire qui affecte la valeur de la série (pandémie, choc boursiers, inflation …)

Pour la suite, nous supposons que la décomposition est additive.

decomposed_data <- decompose(data_ts_train, type="additive")
plot(decomposed_data$trend)

plot(decomposed_data$seasonal)

plot(decomposed_data$random)

boxplot(data_ts ~ cycle(data_ts))

Le choix d’une décomposition additive semble visuellement correct.

Analyse des résidus

L’ACF est le test de l’auto-corrélation par le décalage d’elle même, c’est un test important car il montre si les états précédents (observations décalées) de la série chronologique ont une influence sur l’état actuel.

Si l’auto-corrélation croise la ligne bleue, cela implique qu’un décalage est significativement corrélé avec la série actuelle.

On peut également regarder la distribution des résidus qui ici suit une courbe gaussienne centrée en 0 avec quelques outliers.

checkresiduals(remainder(decomposed_data))
## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.

Représentation de la saisonnalité

En superposant les différentes années on peut observer une saisonnalité qui se répète chaque année.

Ce phénomène est encore plus remarquable quand on supprime la tendance.

On peut noter plusieurs “pic haut” :

  • Décembre / Janvier

  • Mars

  • Juin

  • Septembre (le plus haut)

Et plusieurs “pics faible” :

  • Avril

  • Juillet - Août (la période estivale)

  • Novembre

ggseasonplot(data_ts)

#on supprime la tendance
data_ts_without_trend = diff(data_ts)
SeasonPlot <-  ggseasonplot(data_ts_without_trend) +
  labs(
    title = "Trafic sans la tendance",
    subtitle = "Visualisation de la saisonnalité",
    x = "Mois",
    y = "Nombre de Voyageurs"
  ) +
  geom_line(size = 1.1, alpha = 0.65) +
  theme_fivethirtyeight() +
  theme(axis.title = element_text()) +
  scale_color_brewer(palette = "Paired") +
  theme(axis.title = element_text(), text = element_text(family = "Rubik"))

SeasonPlot

Modèles Espace-état

Pour commencer nos prédictions des 8 premiers mois de l’année 2019, nous allons essayer quelques modèles basiques :

Pour évaluer notre modèle on regarde :

Au formules associés :

La plus populaire est la MAPE

MAPE(y_pred, y_true)

Concrètement si on a une MAPE de valeur 6% cela signifie la différence moyenne entre la prédiction et la valeur réelle est de 6%.

library(forecast)
mean <- meanf(data_ts_train, h=8)
naivem <- naive(data_ts_train, h=8)
driftm <- rwf(data_ts_train, h=8, drif=T)
snaivem <- snaive(data_ts_train, h=8)
## Warning in plot.window(xlim, ylim, log, ...): "plot.conf" n'est pas un paramètre
## graphique
## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "plot.conf" n'est
## pas un paramètre graphique
## Warning in axis(1, ...): "plot.conf" n'est pas un paramètre graphique
## Warning in axis(2, ...): "plot.conf" n'est pas un paramètre graphique
## Warning in box(...): "plot.conf" n'est pas un paramètre graphique

## Warning in plot.window(xlim, ylim, log, ...): "plot.conf" n'est pas un paramètre
## graphique
## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "plot.conf" n'est
## pas un paramètre graphique
## Warning in axis(1, ...): "plot.conf" n'est pas un paramètre graphique
## Warning in axis(2, ...): "plot.conf" n'est pas un paramètre graphique
## Warning in box(...): "plot.conf" n'est pas un paramètre graphique

## Warning in plot.window(xlim, ylim, log, ...): "plot.conf" n'est pas un paramètre
## graphique
## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "plot.conf" n'est
## pas un paramètre graphique
## Warning in axis(1, ...): "plot.conf" n'est pas un paramètre graphique
## Warning in axis(2, ...): "plot.conf" n'est pas un paramètre graphique
## Warning in box(...): "plot.conf" n'est pas un paramètre graphique

#summary(mean)
checkresiduals(mean)

## 
##  Ljung-Box test
## 
## data:  Residuals from Mean
## Q* = 731.64, df = 18, p-value < 2.2e-16
## 
## Model df: 1.   Total lags used: 19
accuracy(mean, data_ts_test)# MAPE : 22%
##                        ME      RMSE       MAE       MPE     MAPE     MASE
## Training set 1.941958e-11  65611.93  55535.08 -3.855657 17.01186 2.161770
## Test set     1.037870e+05 110031.92 103787.04 22.486144 22.48614 4.040036
##                   ACF1 Theil's U
## Training set 0.8254447        NA
## Test set     0.0485288  2.517689
summary(naivem)
## 
## Forecast method: Naive method
## 
## Model Information:
## Call: naive(y = data_ts_train, h = 8) 
## 
## Residual sd: 36679.9508 
## 
## Error measures:
##                    ME     RMSE      MAE         MPE     MAPE     MASE
## Training set 1896.811 36679.95 29013.27 -0.02007386 8.597313 1.129377
##                    ACF1
## Training set -0.2744236
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2019         426097 379089.8 473104.2 354205.6 497988.4
## Feb 2019         426097 359618.7 492575.3 324427.2 527766.8
## Mar 2019         426097 344678.1 507515.9 301577.5 550616.5
## Apr 2019         426097 332082.5 520111.5 282314.2 569879.8
## May 2019         426097 320985.6 531208.4 265343.0 586851.0
## Jun 2019         426097 310953.2 541240.8 249999.8 602194.2
## Jul 2019         426097 301727.5 550466.5 235890.3 616303.7
## Aug 2019         426097 293140.4 559053.6 222757.5 629436.5
checkresiduals(naivem)

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 248.52, df = 19, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 19
accuracy(naivem, data_ts_test) # MAPE : 8.5%
##                     ME     RMSE      MAE         MPE     MAPE     MASE
## Training set  1896.811 36679.95 29013.27 -0.02007386 8.597313 1.129377
## Test set     24357.125 43915.19 38328.62  4.72582164 8.499751 1.491988
##                    ACF1 Theil's U
## Training set -0.2744236        NA
## Test set      0.0485288  1.063155
summary(driftm)
## 
## Forecast method: Random walk with drift
## 
## Model Information:
## Call: rwf(y = data_ts_train, h = 8, drift = T) 
## 
## Drift: 1896.8105  (se 3778.1861)
## Residual sd: 36825.2032 
## 
## Error measures:
##                        ME     RMSE      MAE        MPE     MAPE    MASE
## Training set 2.297696e-11 36630.87 28899.04 -0.5861884 8.591266 1.12493
##                    ACF1
## Training set -0.2744236
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2019       427993.8 380800.4 475187.2 355817.7 500169.9
## Feb 2019       429890.6 362798.7 496982.5 327282.4 532498.8
## Mar 2019       431787.4 349190.1 514384.7 305465.7 558109.1
## Apr 2019       433684.2 337818.7 529549.8 287070.6 580297.9
## May 2019       435581.1 327854.7 543307.4 270827.8 600334.3
## Jun 2019       437477.9 318875.0 556080.7 256090.5 618865.2
## Jul 2019       439374.7 310630.0 568119.3 242476.7 636272.6
## Aug 2019       441271.5 302958.0 579585.0 229739.3 652803.7
checkresiduals(driftm)

## 
##  Ljung-Box test
## 
## data:  Residuals from Random walk with drift
## Q* = 248.52, df = 18, p-value < 2.2e-16
## 
## Model df: 1.   Total lags used: 19
accuracy(driftm, data_ts_test)   # MAPE : 7.5%
##                        ME     RMSE      MAE        MPE     MAPE     MASE
## Training set 2.297696e-11 36630.87 28899.04 -0.5861884 8.591266 1.124930
## Test set     1.582148e+04 41314.60 33586.60  2.7843152 7.582963 1.307399
##                     ACF1 Theil's U
## Training set -0.27442358        NA
## Test set      0.06907259  1.007801
summary(snaivem)
## 
## Forecast method: Seasonal naive method
## 
## Model Information:
## Call: snaive(y = data_ts_train, h = 8) 
## 
## Residual sd: 28666.7301 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE     MAPE MASE      ACF1
## Training set 25337.46 28666.73 25689.63 7.101745 7.207375    1 0.2695124
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2019         415292 378554.1 452029.9 359106.2 471477.8
## Feb 2019         423665 386927.1 460402.9 367479.2 479850.8
## Mar 2019         478207 441469.1 514944.9 422021.2 534392.8
## Apr 2019         443548 406810.1 480285.9 387362.2 499733.8
## May 2019         464162 427424.1 500899.9 407976.2 520347.8
## Jun 2019         457944 421206.1 494681.9 401758.2 514129.8
## Jul 2019         440436 403698.1 477173.9 384250.2 496621.8
## Aug 2019         366272 329534.1 403009.9 310086.2 422457.8
checkresiduals(snaivem)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 35.426, df = 19, p-value = 0.0124
## 
## Model df: 0.   Total lags used: 19
accuracy(snaivem, data_ts_test) # MAPE : 3.6%
##                    ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set 25337.46 28666.73 25689.63 7.101745 7.207375 1.0000000  0.2695124
## Test set     14263.38 22148.43 16960.88 3.053421 3.648407 0.6602226 -0.5427745
##              Theil's U
## Training set        NA
## Test set     0.4792835

Etude du Modèle de Buys-Ballot

Modèle

L’approche de BUYS-BALLOT consiste à introduire des variables indicatrices correspondant à chaque saison définit par le cycle d’observation. Pour les données trimestrielles, on intègre 4 variables indicatrices. Et pour les données mensuelles, on intègre 12 variables indicatrices.

Le modèle doit alors être estimé (sans constante) avec ces variables indicatrices.

Soit le modèle de Buys-Ballot tel que \(Xt = Zt + St + \mu t\)

Nous allons estimé dans un premier temp la tendance \(Zt\), puis dans un second temps la saisonnalité \(St\), tandis que les \(\mu t\) ne pourrons pas être estimer puisqu’il s’agit par définition d’accidents.

Préparation des données

Annees=as.numeric(time(data_ts_train))
ts_DataFrame =data.frame(trafic=data_ts_train,X=as.numeric(Annees))

Estimation de la tendance

Regression <- lm(trafic~X,data = ts_DataFrame)

On utilise la droite de régression obtenu pour prédire la tendance des 8 mois de 2019

tendance=predict(Regression)

#les 8 prochains mois
AnneeMoisNumericFutur=seq(max(Annees)+1/12,length=8,by=1/12)  

tendance2=predict(Regression, newdata=data.frame(X=AnneeMoisNumericFutur)) 

Estimation de la saisonnalité

On récupère les résidus issu de l’estimation de la tendance.

ts_DataFrame$trafic_residual <- residuals(Regression)

On définit la saisonnalité avec les mois.

ts_DataFrame$mois <- round(ts_DataFrame$X - trunc(ts_DataFrame$X),digit=4)

On estime la saisonnalité à partir des résidus du modèle comportant uniquement la tendance.

Regression2 =lm(trafic_residual~0+as.factor(mois),data=ts_DataFrame)

Prédiction de la saisonnalité sur le jeu d’entraînement

prediction2 =predict(Regression2)

Prédiction sur les 8 mois de 2019

MoisNumeric= round(AnneeMoisNumericFutur - trunc(AnneeMoisNumericFutur
                     ),4)
Prediction3 =predict( Regression2, newdata= data.frame(mois=MoisNumeric))

Calculons une région de confiance avec l’erreur d’ajustement

ResidusRegression2=residuals(Regression2)
hist(ResidusRegression2)

1.96*sqrt(var(ResidusRegression2))
## [1] 19226.16

Auto corrélation de la série temporelle

L’auto corrélation de notre série temporelle correspond à la corrélation entre une mesure du trafic \(t\) et les mesures précédentes \(t - k\) ou les mesures suivantes \(t + k\).

L’auto covariance d’une variable \(Xt\) de moyenne \(\mu\) et d’écart type \(\sigma\) à un décalage \(k\) est donné par la formule

\(\gamma_k= E((X_t-\mu)(X_{t+k}-\mu))\)

On en déduit l’auto-corrélation correspondante :

\(\rho_k=\frac{\gamma_k}{\sigma^2}\)

Affichons les auto-corrélations de la séries prédite grâce à un corrélogramme

ACF_Sur_Valeurs_Predites <- acf(prediction2)

Il est normal que la série soit auto corrélé totalement à elle avec un décalage nulle.

On observe une corrélation forte (0.87) avec un décalage (lag) de 12, cela correspond bien à une saisonnalité annuelle.

print(data.frame(ACF_Sur_Valeurs_Predites$lag,ACF_Sur_Valeurs_Predites$acf)[1:13,])
##    ACF_Sur_Valeurs_Predites.lag ACF_Sur_Valeurs_Predites.acf
## 1                             0                   1.00000000
## 2                             1                   0.13556225
## 3                             2                  -0.30353869
## 4                             3                   0.19244891
## 5                             4                  -0.11749420
## 6                             5                  -0.36058135
## 7                             6                  -0.03501791
## 8                             7                  -0.32578902
## 9                             8                  -0.14275858
## 10                            9                   0.15718886
## 11                           10                  -0.25819898
## 12                           11                   0.12067802
## 13                           12                   0.87497657

Recalculons la valeur d’auto-corrélation obtenu en appliquant la formule.

Observons l’application de la formule, en choisissant un décalage de 12

#Constantes
Nombre_Observations=96
decalage=12

#Estimations
moyenneMu=mean(prediction2)
sdSigma=sd(prediction2)


Serie1=prediction2[(decalage+1): 96   ]
Serie2=prediction2[   1 :(96-decalage)]

GammaDecalage12=mean((Serie1-moyenneMu)*(Serie2-moyenneMu))*((Nombre_Observations-decalage)/(Nombre_Observations))

RhoDecalage12=GammaDecalage12/(sdSigma^2)
RhoDecalage12
## [1] 0.8658622

Le résultat obtenu est correct. L’auto corrélation avec un décalage de 12 est donc très forte.

De plus cette auto corrélation étant positive, cela indique une tendance croissante.

la deuxième plus forte corrélation est observé avec un décalage de 5, observons cela graphiquement

plot  ( 1:length(prediction2),   prediction2,type="l")
points((1:length(prediction2))-5,prediction2,type="l",col="red")

Cette corrélation est peu pertinente.

Après avoir étudier les auto-corrélations sur l’ensemble du modèle, Observons les auto-corrélations sur les résidus du modèle de Buys-Ballot.

Une auto corrélation de ces résidus signifierait que les accidents seraient prévisibles, cela n’aurait aucun sens et invaliderait notre modèle.

plot(acf(ResidusRegression2))

Pour notre modèle, il n’y a aucune auto-corrélation significative des résidus (symbolisé par la ligne bleu).

Notre modèle est donc viable.

Comparaison des prédictions et des valeurs réelles

Affichage de la tendance

Affichage du modèle de Buys Ballot

Affichage de la prédiction sur les 8 mois de 2020

Préparation DataFrame pour affichage ggplot

Reproduisons les graphiques avec ggplot2 pour un résultat plus professsionnel.

plotBuysBallot <- Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$BuysBalotModele)
plotBuysBallot

Sauvegarde de l’image pour utilisation ultérieur dans l’application web.

Nous avons réussi à Créer un modèle de prédiction correct avec Buys-Ballot. On remarque que la prédiction semble bien correspondre à la réalité si on fait abstraction du dernier mois où le nombre de voyageurs a bien plus chuté que la prédiction du modèle de Buys-Balot.

Comparons avec un ajustement local réalisé par lissage moyennes mobiles.

Comparaison avec les valeurs observées

Lissage moyenne mobile

Définition

Une moyenne mobile est un filtre linéaire. Il permet de transformer une série chronologique avec comme but d’annuler une composante (tendance ou saison) pour en laisser les autres invariantes tout en réduisant le bruit.

Une moyenne mobile en t est définit comme une combinaison linéaire finie des valeurs de la série correspondantes à des dates entourant t, c’est donc un lissage de la série.

Une moyenne mobile d’ordre m peut être écrite tel que \(\begin{equation} \hat{T}_{t} = \frac{1}{m} \sum_{j=-k}^k y_{t+j}, \end{equation}\)

Avec \(m=2k+1\)

L’estimation de la tendance + saisonnalité au temps t est obtenue en faisant la moyenne des valeurs de la série chronologique dans les k périodes de t. Les observations qui sont proches dans le temps sont également susceptibles d’être proches en valeur estimé.

Par conséquent, la moyenne élimine une partie du caractère aléatoire des données, laissant une composante tendance + saisonnalité lisse.

Choix Moyenne mobiles

Commençons notre analyse avec une moyenne mobile d’ordre 11, décidé de façon arbitraire

library(forecast)
MM_ordre_11 <-  ma(DataAffichageGGplot$trafic[1:96] , 11,centre=FALSE)

La première valeur (non manquante) de cette colonne est la moyenne des onze premières observations càd de janvier à novembre 2011. La deuxième valeur est la moyenne des valeurs de février 2011 à décembre 2011 et ainsi de suite. De façon à que chaque valeur est la moyenne des observations sur 11 mois centrée sur le mois correspondant.

Il est donc normale qu’il manque 5 estimations pour les premiers & derniers mois puisque les données antérieur à 2011 ne sont pas disponibles tout comme les données futur sont inconnus.

Il existe une méthode qui permet l’estimation des points extrêmes que nous utiliserons pour la suite.

Affichage de la moyenne mobile simple d’ordre 11

Affichage_Prediction(DataAffichageGGplot,c(MM_ordre_11,  rep(NA, 8)) )
## Warning: Removed 98 row(s) containing missing values (geom_path).
## Warning: Removed 104 row(s) containing missing values (geom_path).
## Warning: Removed 18 row(s) containing missing values (geom_path).

On remarque que les données estimées (ici en bleu) représente bien la tendance générale de la série, dans le cas d’une moyenne mobile, plus son ordre est élevé, plus la courbe des estimations sera lisse.

Voici pour exemple l’affichage de différentes moyennes mobiles simples.

## Warning: Removed 98 row(s) containing missing values (geom_path).
## Warning: Removed 104 row(s) containing missing values (geom_path).
## Warning: Removed 10 row(s) containing missing values (geom_path).

## Warning: Removed 98 row(s) containing missing values (geom_path).
## Warning: Removed 104 row(s) containing missing values (geom_path).
## Warning: Removed 12 row(s) containing missing values (geom_path).

## Warning: Removed 98 row(s) containing missing values (geom_path).
## Warning: Removed 104 row(s) containing missing values (geom_path).
## Warning: Removed 22 row(s) containing missing values (geom_path).

Nous avons choisi des moyennes mobiles simples d’un ordre impair pour qu’elles soient symétriques. En effet, dans une moyenne mobile d’ordre m=2k+1, on fait la moyenne de l’observation centrale et des k observations de chaque côté. Mais si m était pair, elle ne serait plus symétrique.

Moyenne mobile d’une Moyenne mobile

Il est possible d’appliquer une moyenne mobile à une moyenne mobile. Cela permet de rendre symétrique une moyenne mobile d’ordre pair.

Appliquons cela par une moyenne mobile d’ordre 2 sur une moyenne mobile d’ordre 4

MM_ordre_4 <-  ma(DataAffichageGGplot$trafic[1:96] , 4,centre=FALSE)
MM_Ordre_2x4 <- ma(DataAffichageGGplot$trafic[1:96], 4, centre=TRUE)

head(MM_ordre_4)
## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1]       NA 256279.5 266477.2 282748.5 281141.8 267147.0
head(MM_Ordre_2x4)
## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1]       NA       NA 261378.4 274612.9 281945.1 274144.4

Une MM 2x4 signifie une MM d’ordre 4 suivi d’une MM d’ordre 2. Concrètement, les 5 premières valeurs de la série chronologique sont 245900, 238000, 263227,277991 et 286691. Les deux premiers termes de la MM d’ordre 4 sera

  1. (245900+238000+263227+277991) /4 = 1 025 118/4 = 256279.5

  2. (238000+263227+277991+286691) /4 = 266477.2

Donc le premier terme de la moyenne mobile d’ordre 2, calculé à partir des 2 premiers termes de la MM d’ordre 4, donnera (256279.5 + 266477.2)/2 = 261378.4

On combinant ces deux moyennes mobiles, nous avons obtenu une moyenne mobile centrée d’ordre 4, les valeurs calculées obtenues sont symétriques, en effet, nous avons appliqué cela aux données :

\(\begin{align*} \hat{T}_{t} &= \frac{1}{2}\Big[ \frac{1}{4} (y_{t-2}+y_{t-1}+y_{t}+y_{t+1}) + \frac{1}{4} (y_{t-1}+y_{t}+y_{t+1}+y_{t+2})\Big] \\ &= \frac{1}{8}y_{t-2}+\frac14y_{t-1} + \frac14y_{t}+\frac14y_{t+1}+\frac18y_{t+2}. \end{align*}\)

Ainsi si l’on combine 2 moyennes mobiles, nous devons appliqués une moyenne mobile pair après une moyenne mobile pair et une moyenne mobile impair après une moyenne mobile impair afin que nos résultats soit symétriques.

Conservation & Annulation

Soit la Moyenne mobile 2x12 \(\hat{T}_{t} = \frac{1}{32}y_{t-6} + \frac{1}{12}y_{t-5} + \frac{1}{12}y_{t-4} + \frac{1}{12}y_{t-3} + \frac{1}{12}y_{t-2} + \frac{1}{12}y_{t-1} +\frac{1}{12}y_{t} + \frac{1}{12}y_{t+1} +\frac{1}{12}y_{t+2} + \frac{1}{12}y_{t+3} + \frac{1}{12}y_{t+4} + \frac{1}{12}y_{t+5} + \frac{1}{32}y_{t+6}.\)

Appliquée au trafic qui a une saisonnalité annuelle, chaque mois de l’année a le même poids, car le premier et le dernier terme s’appliquent au même mois des années consécutives. Ainsi, la variation saisonnière est éliminée.

Un autre choix de Moyenne mobile entraînera une estimation soumis à la saisonnalité.

MM_ordre_2x12 <- ma(DataAffichageGGplot$trafic[1:96], 12)
Affichage_Prediction(DataAffichageGGplot,c(MM_ordre_2x12,  rep(NA, 8)) ) + labs(title = "Moyenne Mobile d'ordre 2x12")
## Warning: Removed 98 row(s) containing missing values (geom_path).
## Warning: Removed 104 row(s) containing missing values (geom_path).
## Warning: Removed 20 row(s) containing missing values (geom_path).

La droite de prédiction obtenu ne présente aucune trace de saisonnalité. Un autre choix de moyenne mobile (hormis 12k \(k \in N\)) aurait présenter des traces de saisonnalités.

Nous avons vu que la combinaison de moyennes mobiles produisait des moyennes mobiles pondérés, la combinaison que nous avons utilisés produisait les poids suivants \(\left[\frac{1}{32},\frac{1}{32},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{32}\right]\)

En effet, une moyenne mobile d’ordre m s’écrit tel que \(\hat{T}_t = \sum_{j=-k}^k a_j y_{t+j}\) Avec \(k=(m-1)/2\) et les poids \(\left[a_{-k},\dots,a_k\right]\) Si leur somme vaut 1 et que ces poids sont symétriques (tq \(a_j = a_{-j}\)) alors la courbe obtenu sera plus lisse.

En conclusion, nous avons pu éliminer les saisonnalités de période 12 et conserver la tendance, notre choix d’ordre de moyenne mobile est correct, on remarque que choisir un ordre arbitrairement grand est inutile, il suffit de comparer les moyennes mobiles d’ordre 15 et 12, celle à 12 annule de façon bien plus performante les saisonnalités.

Lissage exponentielle

L’idée de la méthode de lissage exponentielle est une extension de la méthode naîve (naivm vu précédemment). Le principe étant de fournir une prévision qui dépend de moyenne pondérées d’observations passées, les poids diminuant de façon exponentielle à mesure que les observations vieillissent.

Plus simplement, les observations les plus récentes ont un poids plus important que les observations les plus anciennes.

Le paramètre de lissage est décidé par le paramètre α.

On remarque que pour un lissage simple, la prédiction est égale à la valeur de la dernière observation.

Lissage simple

fcst_se <- ses(data_ts_train, h = 8)
summary(fcst_se)
## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = data_ts_train, h = 8) 
## 
##   Smoothing parameters:
##     alpha = 0.2559 
## 
##   Initial states:
##     l = 258126.0245 
## 
##   sigma:  31480.96
## 
##      AIC     AICc      BIC 
## 2430.727 2430.988 2438.420 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE     MAPE     MASE      ACF1
## Training set 7057.127 31151.31 25752.89 1.326234 7.684321 1.002462 0.0711143
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2019       431488.3 391143.9 471832.8 369786.8 493189.9
## Feb 2019       431488.3 389843.9 473132.7 367798.7 495178.0
## Mar 2019       431488.3 388583.3 474393.3 365870.8 497105.9
## Apr 2019       431488.3 387358.8 475617.9 363998.0 498978.7
## May 2019       431488.3 386167.3 476809.4 362175.7 500800.9
## Jun 2019       431488.3 385006.3 477970.4 360400.2 502576.5
## Jul 2019       431488.3 383873.6 479103.1 358667.9 504308.8
## Aug 2019       431488.3 382767.3 480209.4 356975.9 506000.8
checkresiduals(fcst_se)

## 
##  Ljung-Box test
## 
## data:  Residuals from Simple exponential smoothing
## Q* = 144.66, df = 17, p-value < 2.2e-16
## 
## Model df: 2.   Total lags used: 19
plot(fcst_se)
lines(data_ts_test, col="red")

df_se = as.data.frame(fcst_se)
predict_value_se <- df_se$`Point Forecast`
MAPE(predict_value_se, data_ts_test) # MAPE = 7.95
## [1] 7.945782

Optimisation du modèle

La librairie Forecast permet de trouver automatiquement le meilleur degré de lissage exponentielle fournissant la meilleure prédiction.

fit_ets <- ets(data_ts_train) 
print(summary(fit_ets))
## ETS(A,A,A) 
## 
## Call:
##  ets(y = data_ts_train) 
## 
##   Smoothing parameters:
##     alpha = 0.1568 
##     beta  = 1e-04 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 248267.1099 
##     b = 2163.3982 
##     s = -17928.3 -29535.73 9295.935 11005.81 -57117.85 -7708.17
##            38272.64 14592.34 16899.53 34763.15 -7344.204 -5195.15
## 
##   sigma:  11014.45
## 
##      AIC     AICc      BIC 
## 2241.611 2249.458 2285.205 
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -458.6799 10054.77 7831.554 -0.2623253 2.371375 0.3048527
##                    ACF1
## Training set 0.09626331
checkresiduals(fit_ets)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,A,A)
## Q* = 7.1794, df = 3, p-value = 0.06639
## 
## Model df: 16.   Total lags used: 19
fcst_ets <- forecast(fit_ets, h=8)
plot(fcst_ets)
lines(data_ts_test, col="red")

df_ets = as.data.frame(fcst_ets)
predict_value_ets = df_ets$`Point Forecast`
MAPE(predict_value_ets, data_ts_test) #MAPE = 2.84
## [1] 2.839718

Affichage GGplot

DataAffichageGGplot$ModeleLissageExponentielle <- c(fcst_ets$fitted ,predict_value_ets )

Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$ModeleLissageExponentielle)
## Warning: Removed 98 row(s) containing missing values (geom_path).
## Removed 98 row(s) containing missing values (geom_path).

Modèle ARIMA / SAMIRA Automatique

ARIMA : AutoRegressive Integrated Moving Average

Le modèle ARIMA est une combinaison du modèle ARMA* combiné à une différentiation (le Integrated)

*Le modèle ARMA (Autoregressive Moving Average Model) est définit par l’idée de décrire une séries temporelles en deux polynômes. Le premier étant pour l’auto regréssion et le second pour la moyenne mobile. ARMA(p,q) où :

  • p est l’ordre du polynôme auto régressif,

  • q est l’ordre du polynôme de moyenne mobile L’équation est données par avec :

  • φ = les paramètres du modèle autorégressif,

  • θ = les paramètres du modèle de moyenne mobile.

  • c = une constante,

  • Σ = notation de sommation,

  • ε = termes d’erreur (bruit blanc).

Une différentiation correspond à retirer les tendances

  • tendance linéaire : une différenciation

  • tendance quadratique : deux différenciations

Le modèle SARIMA est une combinaison du modèle ARIMA qui prend en compte la composante saisonnière.

Auto.arima prend en compte les saisonnalités, comme on peut le voir dans le modèle sélectionné : (0,1,1)(0,1,1)[12]

# retourne les meilleurs paramètres 
# d=1 enleve la tendance
# D=1 enleve la saisonnalité 
# => avoir des données stationnaires
# trace : voir les résultats
fit_arima <- auto.arima(data_ts_train, d=1, D=1, stepwise = FALSE, approximation = FALSE, trace=TRUE)
## 
##  ARIMA(0,1,0)(0,1,0)[12]                    : 1846.398
##  ARIMA(0,1,0)(0,1,1)[12]                    : 1833.134
##  ARIMA(0,1,0)(0,1,2)[12]                    : 1835.211
##  ARIMA(0,1,0)(1,1,0)[12]                    : 1833.056
##  ARIMA(0,1,0)(1,1,1)[12]                    : 1835.09
##  ARIMA(0,1,0)(1,1,2)[12]                    : Inf
##  ARIMA(0,1,0)(2,1,0)[12]                    : 1835.207
##  ARIMA(0,1,0)(2,1,1)[12]                    : 1837.012
##  ARIMA(0,1,0)(2,1,2)[12]                    : 1836.461
##  ARIMA(0,1,1)(0,1,0)[12]                    : 1814.951
##  ARIMA(0,1,1)(0,1,1)[12]                    : 1801.155
##  ARIMA(0,1,1)(0,1,2)[12]                    : 1803.362
##  ARIMA(0,1,1)(1,1,0)[12]                    : 1803.592
##  ARIMA(0,1,1)(1,1,1)[12]                    : 1803.361
##  ARIMA(0,1,1)(1,1,2)[12]                    : Inf
##  ARIMA(0,1,1)(2,1,0)[12]                    : 1805.004
##  ARIMA(0,1,1)(2,1,1)[12]                    : 1805.397
##  ARIMA(0,1,1)(2,1,2)[12]                    : Inf
##  ARIMA(0,1,2)(0,1,0)[12]                    : 1816.915
##  ARIMA(0,1,2)(0,1,1)[12]                    : 1803.033
##  ARIMA(0,1,2)(0,1,2)[12]                    : 1805.296
##  ARIMA(0,1,2)(1,1,0)[12]                    : 1805.702
##  ARIMA(0,1,2)(1,1,1)[12]                    : 1805.295
##  ARIMA(0,1,2)(1,1,2)[12]                    : Inf
##  ARIMA(0,1,2)(2,1,0)[12]                    : 1807.026
##  ARIMA(0,1,2)(2,1,1)[12]                    : 1807.441
##  ARIMA(0,1,3)(0,1,0)[12]                    : 1817.787
##  ARIMA(0,1,3)(0,1,1)[12]                    : Inf
##  ARIMA(0,1,3)(0,1,2)[12]                    : Inf
##  ARIMA(0,1,3)(1,1,0)[12]                    : Inf
##  ARIMA(0,1,3)(1,1,1)[12]                    : Inf
##  ARIMA(0,1,3)(2,1,0)[12]                    : Inf
##  ARIMA(0,1,4)(0,1,0)[12]                    : 1820.052
##  ARIMA(0,1,4)(0,1,1)[12]                    : Inf
##  ARIMA(0,1,4)(1,1,0)[12]                    : Inf
##  ARIMA(0,1,5)(0,1,0)[12]                    : Inf
##  ARIMA(1,1,0)(0,1,0)[12]                    : 1825.579
##  ARIMA(1,1,0)(0,1,1)[12]                    : 1812.512
##  ARIMA(1,1,0)(0,1,2)[12]                    : 1814.657
##  ARIMA(1,1,0)(1,1,0)[12]                    : 1813.2
##  ARIMA(1,1,0)(1,1,1)[12]                    : 1814.614
##  ARIMA(1,1,0)(1,1,2)[12]                    : 1816.192
##  ARIMA(1,1,0)(2,1,0)[12]                    : 1815.227
##  ARIMA(1,1,0)(2,1,1)[12]                    : Inf
##  ARIMA(1,1,0)(2,1,2)[12]                    : 1817.796
##  ARIMA(1,1,1)(0,1,0)[12]                    : 1816.841
##  ARIMA(1,1,1)(0,1,1)[12]                    : 1802.853
##  ARIMA(1,1,1)(0,1,2)[12]                    : 1805.117
##  ARIMA(1,1,1)(1,1,0)[12]                    : 1805.653
##  ARIMA(1,1,1)(1,1,1)[12]                    : Inf
##  ARIMA(1,1,1)(1,1,2)[12]                    : Inf
##  ARIMA(1,1,1)(2,1,0)[12]                    : Inf
##  ARIMA(1,1,1)(2,1,1)[12]                    : Inf
##  ARIMA(1,1,2)(0,1,0)[12]                    : 1819.234
##  ARIMA(1,1,2)(0,1,1)[12]                    : 1805.22
##  ARIMA(1,1,2)(0,1,2)[12]                    : 1807.539
##  ARIMA(1,1,2)(1,1,0)[12]                    : 1807.381
##  ARIMA(1,1,2)(1,1,1)[12]                    : 1807.538
##  ARIMA(1,1,2)(2,1,0)[12]                    : 1808.925
##  ARIMA(1,1,3)(0,1,0)[12]                    : 1820.05
##  ARIMA(1,1,3)(0,1,1)[12]                    : 1806.055
##  ARIMA(1,1,3)(1,1,0)[12]                    : 1808.732
##  ARIMA(1,1,4)(0,1,0)[12]                    : Inf
##  ARIMA(2,1,0)(0,1,0)[12]                    : 1824.435
##  ARIMA(2,1,0)(0,1,1)[12]                    : 1811.07
##  ARIMA(2,1,0)(0,1,2)[12]                    : 1813.287
##  ARIMA(2,1,0)(1,1,0)[12]                    : 1811.619
##  ARIMA(2,1,0)(1,1,1)[12]                    : 1813.247
##  ARIMA(2,1,0)(1,1,2)[12]                    : Inf
##  ARIMA(2,1,0)(2,1,0)[12]                    : 1813.821
##  ARIMA(2,1,0)(2,1,1)[12]                    : 1815.872
##  ARIMA(2,1,1)(0,1,0)[12]                    : Inf
##  ARIMA(2,1,1)(0,1,1)[12]                    : Inf
##  ARIMA(2,1,1)(0,1,2)[12]                    : Inf
##  ARIMA(2,1,1)(1,1,0)[12]                    : Inf
##  ARIMA(2,1,1)(1,1,1)[12]                    : Inf
##  ARIMA(2,1,1)(2,1,0)[12]                    : Inf
##  ARIMA(2,1,2)(0,1,0)[12]                    : Inf
##  ARIMA(2,1,2)(0,1,1)[12]                    : Inf
##  ARIMA(2,1,2)(1,1,0)[12]                    : Inf
##  ARIMA(2,1,3)(0,1,0)[12]                    : Inf
##  ARIMA(3,1,0)(0,1,0)[12]                    : 1823.646
##  ARIMA(3,1,0)(0,1,1)[12]                    : 1808.49
##  ARIMA(3,1,0)(0,1,2)[12]                    : 1810.542
##  ARIMA(3,1,0)(1,1,0)[12]                    : 1808.594
##  ARIMA(3,1,0)(1,1,1)[12]                    : 1810.321
##  ARIMA(3,1,0)(2,1,0)[12]                    : 1810.708
##  ARIMA(3,1,1)(0,1,0)[12]                    : Inf
##  ARIMA(3,1,1)(0,1,1)[12]                    : Inf
##  ARIMA(3,1,1)(1,1,0)[12]                    : Inf
##  ARIMA(3,1,2)(0,1,0)[12]                    : Inf
##  ARIMA(4,1,0)(0,1,0)[12]                    : 1823.996
##  ARIMA(4,1,0)(0,1,1)[12]                    : 1810.199
##  ARIMA(4,1,0)(1,1,0)[12]                    : 1810.845
##  ARIMA(4,1,1)(0,1,0)[12]                    : Inf
##  ARIMA(5,1,0)(0,1,0)[12]                    : 1825.055
## 
## 
## 
##  Best model: ARIMA(0,1,1)(0,1,1)[12]
print(summary(fit_arima))
## Series: data_ts_train 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.7675  -0.5465
## s.e.   0.0977   0.1295
## 
## sigma^2 = 138827719:  log likelihood = -897.43
## AIC=1800.85   AICc=1801.16   BIC=1808.11
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 805.2378 10822.93 7747.841 0.2158443 2.213481 0.3015941 0.03453594
checkresiduals(fit_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 10.898, df = 17, p-value = 0.8618
## 
## Model df: 2.   Total lags used: 19
fcst_arima <- forecast(fit_arima, h=8)
plot(fcst_arima)
lines(data_ts_test, col='red')

df_arima = as.data.frame(fcst_arima)
predict_value_arima = df_arima$`Point Forecast`
MAPE(predict_value_arima, data_ts_test)
## [1] 2.72507

Affichage GGplot

DataAffichageGGplot$ModeleArima <-  c(fit_arima$fitted ,predict_value_arima )
Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$ModeleArima)
## Warning: Removed 98 row(s) containing missing values (geom_path).
## Removed 98 row(s) containing missing values (geom_path).

Prophet

Prophet est un outil de prévision de série chronologique open-source mis en production par Facebook.

L’idée est décomposer la série temporelle sous la forme :

y(t) = g(t) + s(t) + h(t) + e(t)

avec respectivement :

  • g(t) : la tendance (linéaire ou logistique)
  • s(s) : une ou plusieurs composantes saisonnières (annuelle, hebdomadaire ou quotidienne)
  • h(t) : l’effet des vacances ou de jours spécifiques qui pourront être paramétrés
  • e(t) : l’erreur, bruit aléatoire qui mesure l’écart entre le modèle et les données réelles

Ce modèle est particulièrement utile quand on a de forte saisonnalités et permet de prendre en compte les vacances (et ses effets) par pays.

Préparation des données

  • format des dates

  • colonne du temps name : ds

  • colonne de la mesure quantitative name : y

library(prophet)
## Le chargement a nécessité le package : Rcpp
## Le chargement a nécessité le package : rlang
## 
## Attachement du package : 'rlang'
## Les objets suivants sont masqués depuis 'package:purrr':
## 
##     %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
##     flatten_lgl, flatten_raw, invoke, splice
library(zoo)
## 
## Attachement du package : 'zoo'
## Les objets suivants sont masqués depuis 'package:base':
## 
##     as.Date, as.Date.numeric
data_train$ds <- as.Date( as.yearmon(time(data_ts_train)))

On créer notre modèle et on prédit sur les 8 mois de 2019.

On affiche les courbes.

On regarde la valeur de la MAPE : 3.2.

model_prophet <- prophet(data_train)
## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
forecast_prophet <- make_future_dataframe(model_prophet, periods = 8, freq = 'month')
AAPLfc <- predict(model_prophet, forecast_prophet)
tail(AAPLfc[c("ds", "yhat", "yhat_lower", "yhat")])
##             ds     yhat yhat_lower   yhat.1
## 99  2019-03-01 493929.1   482512.9 493929.1
## 100 2019-04-01 477074.4   465193.4 477074.4
## 101 2019-05-01 475880.5   462887.3 475880.5
## 102 2019-06-01 503278.4   491122.3 503278.4
## 103 2019-07-01 457437.7   445628.0 457437.7
## 104 2019-08-01 412673.2   401309.4 412673.2
dyplot.prophet(model_prophet, AAPLfc)
## Warning: `select_()` was deprecated in dplyr 0.7.0.
## Please use `select()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
data_pp <- subset(AAPLfc, select=c("yhat"))
data_pp_ts <- ts(data_pp, start=2011, frequency=12)
data_pp_ts_w <- window(data_pp_ts, start= c(2019,1), end = c(2019,8))
MAPE(data_pp_ts_w, data_ts_test) #3.2
## [1] 3.189021

Affichage GGplot

DataAffichageGGplot$ProphetModele <- AAPLfc$yhat
Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$ProphetModele)
## Warning: Removed 98 row(s) containing missing values (geom_path).
## Removed 98 row(s) containing missing values (geom_path).

LSTM

Un réseau LSTM (Long Short Term Memory) est un type particulier de réseau de neurones récurrent qui est capable d’apprendre les dépendances à long terme. Il a été introduit pas Hochreiter & Schmidhuber en 1997.

L’idée d’un réseaux récurrent est de mettre en place des boucles qui permettent aux informations de persister. Ils sont représentés par une chaîne de modules répétitifs de réseau de neurones.

Chaque module est traversé un à un en faisant évoluer son état qui peut se voir ajouter ou supprimer de l’information, ces changements sont régulés par ce qu’on appelle des “portes”.

Elles sont définit par une couche de neurones sigmoïde et d’opération de multiplication ponctuelle qui laisse passer ou non l’information. La fonction Sigmoïde produit des nombres entre 1 et 0 ce qui définit la quantité de chaque composant qui doit passer.

Préparation des données

On veut prendre l’année précédente pour apprendre , on fixe donc lag à 12, en réalité cela correspond à faire 12 - 1.

Puis on transforme en array 3D car le modèle LSTM prend un tensor de format 3D [samples, timesteps, features]

Construction de notre réseaux de neuronnes

lstm_model <- keras_model_sequential()
## Loaded Tensorflow version 2.8.0
lstm_model %>%
  layer_lstm(units = 50, # taille du layer
       batch_input_shape = c(1, 12, 1), # batch size, timesteps, features
       return_sequences = TRUE,
       stateful = TRUE) %>%
  # pourcentage de neuronnes qu'on détruit à chaque itération 
  layer_dropout(rate = 0.5) %>%
  layer_lstm(units = 50,
        return_sequences = TRUE,
        stateful = TRUE) %>%
  layer_dropout(rate = 0.5) %>%
  time_distributed(keras::layer_dense(units = 1))

lstm_model %>%
    compile(loss = 'mae', optimizer = 'adam', metrics = 'accuracy')

summary(lstm_model)
## Model: "sequential"
## ________________________________________________________________________________
##  Layer (type)                       Output Shape                    Param #     
## ================================================================================
##  lstm_1 (LSTM)                      (1, 12, 50)                     10400       
##  dropout_1 (Dropout)                (1, 12, 50)                     0           
##  lstm (LSTM)                        (1, 12, 50)                     20200       
##  dropout (Dropout)                  (1, 12, 50)                     0           
##  time_distributed (TimeDistributed)  (1, 12, 1)                     51          
## ================================================================================
## Total params: 30,651
## Trainable params: 30,651
## Non-trainable params: 0
## ________________________________________________________________________________

Apprentissage sur le jeu d’entraînement

lstm_model %>% fit(
    x = x_train_arr,
    y = y_train_arr,
    batch_size = 1,
    epochs = 20,
    verbose = 0,
    shuffle = FALSE #pour preserver nos séquences
)

Prédiction sur le jeu de Test

On a 50 résultats / prédictions par input donc on transforme pour une seule prédiction

lstm_forecast <- lstm_model %>%
    predict(x_pred_arr, batch_size = 1) %>%
    .[, , 1]
 
# rescale en format basique
lstm_forecast <- lstm_forecast * scale_factors[2] + scale_factors[1]
lstm_forecast

fitted <- predict(lstm_model, x_train_arr, batch_size = 1) %>%
     .[, , 1]

if (dim(fitted)[2] > 1) {
  
    fit <- c(fitted[, 1], fitted[dim(fitted)[1], 2:dim(fitted)[2]])
} else {
    fit <- fitted[, 1]
}

# rescale final de nos données
fitted <- fit * scale_factors[2] + scale_factors[1]
fitted
fitted <- c(rep(NA, lag), fitted)
fitted
length(fitted)

Affichage et évaluation

Nos résultats ne sont pas très convaincant même en faisant varier nos paramètres, ceci peut être explicable par le manque d’observations.

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2019 469206.2 471857.1 437201.1 446503.8 458682.0 457769.1 478484.8 472561.8
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug
## 2019 443700 441499 480649 463680 453372 505190 445332 370211

Affichage GGplot

DataAffichageGGplot$LSTM_Modele <- fitted
Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$LSTM_Modele)
## Warning: Removed 98 row(s) containing missing values (geom_path).
## Removed 98 row(s) containing missing values (geom_path).
## Warning: Removed 12 row(s) containing missing values (geom_path).

Comparaison Modèle

Choix des indicateurs

Nous devons pouvoir comparer nos modèles entre eux, nous souhaitons pouvoir observer leurs performances globale, ainsi que spécifiquement leurs performances sur le jeu d’entraînement et le jeu de test (les 8 mois de 2019 à prédire).

Nous utiliserons une mesure largement utilisé lors de nos précédents travaux : le \(R^2 \in [0 \; ; 1]\) , cette mesure indique la proportion de la variance expliquée par le modèle.

  • 0 % le modèle n’explique par la variable Y

  • 100 % le modèle explique la variabilité de Y lié à la liaison linéaire des variables explicatives entièrement

Rappelons son explication mathématiques , soit SCE la somme des distances au carré entre chaque valeur prédite par le modèle \({\widehat y_i}\) et la moyenne des réponses \(\overline{y}\)

\(\text{SCE} = \sum_{i=1}^{N}(\hat{y_i} – \overline{y})^2\)

Nous obtenons alors la part de dispersion expliquée par le modèle.

Puis, nous calculons la dispersion totale des données nommé SCT

\(\text{SCT} = \sum_{i=1}^{N}(y_{i } – \overline{y})^2\)

Avec \(y_i\) une valeur prise par une variable expliquée

Nous obtenons alors le R² par la combinaisons des calculs précédents

\(R^2 = \frac{SCE}{SCT}\)

Pour compléter cette mesure, nous utiliserons donc l’erreur absolue moyenne en pourcentage (MAPE en anglais) Il s’agit de la moyenne des écarts en valeur absolue par rapport aux valeurs observées.

C’est donc un pourcentage et par conséquent un indicateur pratique de comparaison.

Exemple de mise en production

Imaginons que l’on utilise nos modèles pour prédire le nombre de voyageurs, si l’on prévoit trop de passagers alors les moyens mis en place pour les accueillir sont partiellement utilisés, il en résulte un coût de 5 euros par passagers.

Si l’on a prédit moins de passagers que la réalité, notre compagnie doit commander en livraison urgente, il en résulte un coût de 10 euros par passagers.

Observons cette mise en production sur nos modèles.

Calculs des performances

Dans les concerts, chaque spectateur chante comme une casserole, pourtant la foule chante systématiquement juste, suivant cet intuition réalisons une prédiction qui est une moyenne des différents modèles réalisés.

DataAffichageGGplot$MoyenneDesModele <- rowMeans(DataAffichageGGplot[,5:9])

Performances Modèles

noms_modeles <- c("Buys Ballot ","Lissage exponentielle ","Arima","Prophet","LSTM","Moyenne Des Modèles")

for (x in seq_along(noms_modeles)){
  print("----------------------------")
  print(noms_modeles[x])
  
  predictionglobale <- DataAffichageGGplot[13:104,4+x]
  predictionEntrainement <- DataAffichageGGplot[13:96,4+x]
  predictionTest <- DataAffichageGGplot[97:104,4+x]
    
  
  cat("MAPE globale : ",MAPE(DataAffichageGGplot$trafic[13:104], predictionglobale),"\n")
cat("R carré globale : ",Rcarre(DataAffichageGGplot$trafic[13:104], predictionglobale),"\n")

cat("MAPE Entrainement : ",MAPE(DataAffichageGGplot$trafic[13:96], predictionEntrainement),"\n")
cat("R carré Entrainement : ",Rcarre(DataAffichageGGplot$trafic[13:96], predictionEntrainement),"\n")

cat("MAPE Test : ",MAPE(DataAffichageGGplot$trafic[97:104], predictionTest),"\n")
cat("R carré Test : ",Rcarre(DataAffichageGGplot$trafic[97:104], predictionTest),"\n")

MiseEnProductionDuModele <-  CoutDesErreurs(DataAffichageGGplot$trafic[97:104], predictionTest)
cat("Nombre de Passagers prévus en plus : ",MiseEnProductionDuModele[1],"\nNombre de Passagers prévus en moins : ",MiseEnProductionDuModele[2],"\nCout des erreurs en Euros : ",MiseEnProductionDuModele[3] )


cat("\nDifférence en nbre de voyageurs sur le dernier mois : ",predictionglobale[92] - DataAffichageGGplot$trafic[104],"\n")
  
}
## [1] "----------------------------"
## [1] "Buys Ballot "
## MAPE globale :  2.205698 
## R carré globale :  0.9740106 
## MAPE Entrainement :  2.142334 
## R carré Entrainement :  0.9747808 
## MAPE Test :  2.871024 
## R carré Test :  0.8036533 
## Nombre de Passagers prévus en plus :  4848.58 
## Nombre de Passagers prévus en moins :  -91908.62 
## Cout des erreurs en Euros :  1016058
## Différence en nbre de voyageurs sur le dernier mois :  36037.92 
## [1] "----------------------------"
## [1] "Lissage exponentielle "
## MAPE globale :  2.270427 
## R carré globale :  0.9725982 
## MAPE Entrainement :  2.200386 
## R carré Entrainement :  0.9736448 
## MAPE Test :  3.005848 
## R carré Test :  0.7863116 
## Nombre de Passagers prévus en plus :  5276.056 
## Nombre de Passagers prévus en moins :  -96021.53 
## Cout des erreurs en Euros :  1065736
## Différence en nbre de voyageurs sur le dernier mois :  38627.88 
## [1] "----------------------------"
## [1] "Arima"
## MAPE globale :  2.551599 
## R carré globale :  0.9655633 
## MAPE Entrainement :  2.526595 
## R carré Entrainement :  0.9629896 
## MAPE Test :  2.814135 
## R carré Test :  0.8420685 
## Nombre de Passagers prévus en plus :  17397.14 
## Nombre de Passagers prévus en moins :  -81186.49 
## Cout des erreurs en Euros :  1159808
## Différence en nbre de voyageurs sur le dernier mois :  24913.57 
## [1] "----------------------------"
## [1] "Prophet"
## MAPE globale :  2.0477 
## R carré globale :  0.9744149 
## MAPE Entrainement :  1.919235 
## R carré Entrainement :  0.9777721 
## MAPE Test :  3.396584 
## R carré Test :  0.7327864 
## Nombre de Passagers prévus en plus :  1911.595 
## Nombre de Passagers prévus en moins :  -112557.7 
## Cout des erreurs en Euros :  1163809
## Différence en nbre de voyageurs sur le dernier mois :  42462.18 
## [1] "----------------------------"
## [1] "LSTM"
## MAPE globale :  3.439605 
## R carré globale :  0.944196 
## MAPE Entrainement :  3.5023 
## R carré Entrainement :  0.9383416 
## MAPE Test :  2.781304 
## R carré Test :  0.7919626 
## Nombre de Passagers prévus en plus :  48365.29 
## Nombre de Passagers prévus en moins :  -51742.04 
## Cout des erreurs en Euros :  1484726
## Différence en nbre de voyageurs sur le dernier mois :  19874.38 
## [1] "----------------------------"
## [1] "Moyenne Des Modèles"
## MAPE globale :  2.210271 
## R carré globale :  0.9741124 
## MAPE Entrainement :  2.151987 
## R carré Entrainement :  0.9742663 
## MAPE Test :  2.822246 
## R carré Test :  0.8218652 
## Nombre de Passagers prévus en plus :  12776.57 
## Nombre de Passagers prévus en moins :  -83900.11 
## Cout des erreurs en Euros :  1094533
## Différence en nbre de voyageurs sur le dernier mois :  32383.19

Sauvegardons le modèle de Buys Ballot pour une utilisation ultérieur dans une application Web

saveRDS(Regression, file = "./Shiny/RegressionAnnees.rda")
saveRDS(Regression2, file = "./Shiny/RegressionMois.rda")

Conclusion